/*
*		
*		A compilacao e feita por: gcc -o gauss_pthread gauss_pthread.c -lgmp -lm -lpthread
*
*/


#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <pthread.h> 
#include <time.h>
#include <gmp.h> 

#define N 25
#define BASE 10
#define MANTISSA 1024*1024*32

#define N_THREADS 2* N


// define as funcoes
void init_mpf_b(mpf_t rop);
void media_simples(mpf_t rop, mpf_t a, mpf_t b);
void media_geometrica(mpf_t rop, mpf_t a, mpf_t b); 
void equacao_t(mpf_t t, mpf_t x, mpf_t y, mpf_t a);
void equacao_pi(mpf_t pi, mpf_t a, mpf_t b, mpf_t t);

//cria os semaforos
pthread_mutex_t MUTEX;
pthread_cond_t isItOk;

// cria as variaveis
mpf_t a;
mpf_t b;
mpf_t t;
mpf_t x;
mpf_t y;
mpf_t pi;


// thread A recebe o id da thread que foi criada 
// calcula a media simples, a geometrica
// define os valores de a, b e y
void *ThreadA(void *tid)
{
    int rc;

    pthread_mutex_lock(&MUTEX);
	
    //a = (a + b)/2.0;
    media_simples(a, a, b);     

    pthread_mutex_unlock(&MUTEX);
    pthread_cond_signal(&isItOk);  
   
    pthread_cond_wait(&isItOk, &MUTEX);
	
	//b = sqrt(b*y);
    media_geometrica(b, b, y);  
	//y = a;
    mpf_set(y, a);             

    pthread_mutex_unlock(&MUTEX);
    
    pthread_exit(NULL);
}


// thread B recebe o id da thread que foi criada 
// calcula a equacao de t
// define os valores de t e x
void *ThreadB(void *tid)
{
	//t = t -   x * pow((y-a),2); 
	equacao_t(t, x, y, a);      
    //x *= 2; 
	mpf_mul_ui(x, x, 2);        
   
    pthread_cond_signal(&isItOk);  
    pthread_mutex_unlock(&MUTEX);
    
	pthread_exit(NULL);
}


int main(int argc, char* argv[])
{
    unsigned long i, j = 0, p = 0, q = 0; // variaveis auxiliares

    pthread_t threads[N_THREADS];

    pthread_attr_t attr;
    void *exit_status;

    time_t t_zero, t_final;
    double delta_t;

    // Marca o tempo
    time(&t_zero);


    // Colocar Pthreads como Joinable    
    pthread_attr_init(&attr);
    pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
    
	// inicializar Mutexes
    pthread_cond_init(&isItOk, NULL);
    pthread_mutex_init(&MUTEX, NULL);

    // Float com precisao arbitraria
    mpf_set_default_prec(MANTISSA);

	
	//inicializa as variaveis
    mpf_init(a);
    mpf_init(b);
    mpf_init(t);
    mpf_init(x);
    mpf_init(y);
    mpf_init(pi);

    //	a = 1.0;
	//	b = (1.0/(sqrt(2)));
	//	t = (1.0/4.0);
	//	x = 1.0;
	//	y = 1.0; // Independente
    mpf_set_d(t, (1.0/4.0));
    mpf_set_ui(a, 1);
    mpf_set_ui(y, 1);
    mpf_set_ui(x, 1);
    init_mpf_b(b);

	
    /*
     *  A(n+1) = A(n)+B(n)/2
     *  B(n+1) = sqrt(A(n)+B(n))
     *  T(n+1) = T(n) - P(n)(A(n) - A(n+1))^2
     *  P(n+1) = 2*P(n)
     *  pi = ((A(n)+B(n))^2)/4T(n)
     * 
     * */

    for(i = 0; i < N; i++)
	{
        p = j;
		// cria thread A
        pthread_create(&threads[j], &attr, ThreadA, (void*) q);    
        j++; q++;

		//espera thread A acabar o calculo
        pthread_cond_wait(&isItOk, &MUTEX);
		
		// cria thread B
        pthread_create(&threads[j], &attr, ThreadB, (void*) q);
        j++; q++;
        
		// sincroniza threads
        for(j = p;j < q; j++)
		   {
             pthread_join(threads[j], NULL);
           }
    }
    
	
	//pi = pow((a+b),2)/(4.0*t);
    equacao_pi(pi, a, b, t);    
    
    // Libera a memoria
    mpf_clear(a);
    mpf_clear(b);
    mpf_clear(y);
    mpf_clear(t);

	
	//mostra o valor de pi
    mpf_out_str(NULL, BASE, 0, pi);
	
	// Libera a memoria
    mpf_clear(pi);
    
	// destroi os semaforos
    pthread_attr_destroy(&attr);
    pthread_mutex_destroy(&MUTEX);
    pthread_cond_destroy(&isItOk);

	//mostra o tempo
    time(&t_final);
    delta_t = difftime(t_final, t_zero);
    printf("\nTempo = %.2lf segundos\n", delta_t);

	// espera todas as threads acabarem de rodar antes da main terminar
    pthread_exit(NULL);   
    
    return 0;
}


// Inicializa o valor de b, onde a variavel b e passada por parametro
// b = 1/sqrt(2)
void init_mpf_b(mpf_t rop)
{
    mpf_set_default_prec(MANTISSA);
	mpf_t x, y;

	mpf_init(x);
	mpf_init(y);
	
	// x = 1
	mpf_set_ui(x, 1);
	// y = sqtr(2)
	mpf_sqrt_ui(y, 2);
	// x = 1/sqrt(2)
	mpf_div(x, x, y);

	// return x
	mpf_set(rop, x);

	//libera memoria
	mpf_clear(x);
	mpf_clear(y);
}


// calcula a media simples entre a e b, onde as variaveis a,b e 'rop' sao passadas por parametro
// (a + b) / 2
void media_simples(mpf_t rop, mpf_t a, mpf_t b)
{
	mpf_t aux;	
	mpf_init(aux);

	//aux = a + b;
	mpf_add(aux, a, b);
	
	//rop = aux/2
	mpf_div_ui(rop, aux, 2);

	//libera memoria
	mpf_clear(aux);
}


// calcula a media simples entre a e b, onde as variaveis a,b e 'rop' sao passadas por parametro
// sqrt(a + b)
void media_geometrica(mpf_t rop, mpf_t a, mpf_t b)
{	
	mpf_t aux;
	mpf_init(aux);

	//aux = a + b;
	mpf_mul(aux, a, b);
	
	//rop = sqrt(aux)
	mpf_sqrt(rop, aux);

	//libera memoria
	mpf_clear(aux);
}


// calcula o valor de t onde as variaveis t, x, y e a sao passadas por parametro
// t = t - x * (y-a)^2;
void equacao_t(mpf_t t, mpf_t x, mpf_t y, mpf_t a)
{
	mpf_t aux;
	mpf_init(aux);
	

	// aux = y - a
	mpf_sub(aux, y, a);
	// aux = (aux)^2
	mpf_pow_ui(aux, aux, 2);
	// aux = aux * x
	mpf_mul(aux, aux, x);
	//aux = t - aux
	mpf_sub(t, t, aux);

	//libera memoria
	mpf_clear(aux);
}


// calcula o valor de pi onde as variaveis pi, a, b e t sao passadas por parametro
// pi = ((a + b)^2)/(4.0*t)
void equacao_pi(mpf_t pi, mpf_t a, mpf_t b, mpf_t t)
{
	mpf_t aux;
	mpf_init(aux);
		
	// aux = a + b
	mpf_add(aux, a, b);
	// aux = aux^2
	mpf_pow_ui(aux, aux, 2);
	// aux = aux / t
	mpf_div(aux, aux, t);
	// pi = aux/4
	mpf_div_ui(pi, aux, 4);

	//libera memoria
	mpf_clear(aux);
}
